#include "TMath.h"

const Double_t EPSILON = 1e-12;

Double_t conv(const Double_t x, const Double_t amp, const Double_t mean, const Double_t sigma, const Double_t lambda)
{
  const Double_t x0 = x-mean;
  Double_t sum=0, sumW=0;

  for(Double_t nlambda=0; nlambda<5; nlambda+=0.5){
    const Double_t weight = TMath::Exp(-nlambda);
    sum += TMath::Gaus( x0 - nlambda*lambda, 0, sigma, kTRUE)*weight;
    sumW+=weight;
  }
  
  if (sumW>0) return amp*sum/sumW;
  return 0;
}

Double_t conv1(const Double_t xx, const Double_t rawamp, const Double_t rawmean, const Double_t rawsigma, const Double_t rawlambda)
{
  //true
  /*
  const Int_t nsample = 20000;
  const Double_t mean = 100;
  const Double_t sigma = 7;
  const Double_t lambda = 10.;
  */
  //step = 0.1
  /*
 FCN=138.085 FROM MIGRAD    STATUS=CONVERGED     188 CALLS         189 TOTAL
                     EDM=1.60535e-10    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   0.7 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           1.98621e+04   1.41066e+02  -7.70558e-02  -1.28797e-07
   2  p1           1.00391e+02   1.15203e-01  -1.77436e-04  -1.79985e-04
   3  p2           6.97207e+00   8.40582e-02  -6.13961e-05   3.19997e-05
   4  p3           1.02199e+01   1.35197e-01   1.36232e-04  -1.79411e-04
   */
  //step=0.01
  /*
 FCN=142.592 FROM MIGRAD    STATUS=CONVERGED     223 CALLS         224 TOTAL
                     EDM=2.30397e-07    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   2.5 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           1.98576e+04   1.40569e+02  -1.37550e+00   2.05244e-06
   2  p1           9.99193e+01   1.19409e-01   2.21809e-04   1.14837e-02
   3  p2           6.95751e+00   8.47136e-02  -6.09743e-06  -9.53140e-03
   4  p3           1.02606e+01   1.36194e-01   7.25162e-05   2.96274e-03
   */
  //step=1
  /*
 FCN=141.1 FROM MIGRAD    STATUS=CONVERGED     231 CALLS         232 TOTAL
                     EDM=2.68669e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           1.98593e+04   1.40928e+02   8.20044e-01  -1.13907e-06
   2  p1           1.04131e+02   8.50050e-02   4.22977e-04   3.18552e-04
   3  p2           7.49839e+00   7.25817e-02   3.61602e-04   2.31573e-05
   4  p3           1.03118e+01   1.36127e-01   6.81078e-04   5.34253e-03
   */
  //step=0.5
  /*
 FCN=142.96 FROM MIGRAD    STATUS=CONVERGED     171 CALLS         172 TOTAL
                     EDM=4.73326e-08    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   1.5 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           1.98572e+04   1.39789e+02  -7.15466e-01  -2.42201e-07
   2  p1           1.02217e+02   9.66143e-02  -3.23427e-04  -2.24813e-03
   3  p2           7.10800e+00   8.18227e-02  -5.50021e-04   2.68938e-03
   4  p3           1.02684e+01   1.33740e-01   1.17429e-03  -2.85964e-03
   */

  //maxn=20, has sumW, step 0.1
  /*
 FCN=169.765 FROM MIGRAD    STATUS=CONVERGED     187 CALLS         188 TOTAL
                     EDM=8.87966e-10    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   1.0 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           1.99844e+05   4.51638e+02  -1.71982e-01  -1.19477e-07
   2  p1           1.00592e+02   3.60566e-02   2.20781e-05  -7.82818e-04
   3  p2           7.05192e+00   2.60868e-02  -2.77496e-06   9.48287e-06
   4  p3           9.89452e+00   4.04602e-02  -9.20382e-05  -7.02322e-04
  */
  //maxn = 20, use mean, step 0.1
  /*
    FCN=169.765 FROM MIGRAD    STATUS=CONVERGED     358 CALLS         359 TOTAL
                     EDM=1.76185e-08    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   3.8 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           1.88170e+05   8.82640e+02  -1.10782e+00  -5.11834e-07
   2  p1           1.00592e+02   3.60611e-02   5.86955e-05   3.86697e-03
   3  p2           7.05192e+00   2.61887e-02   3.38056e-05  -3.69544e-03
   4  p3           9.89452e+00   4.04989e-02  -7.54308e-05   8.85499e-03
   */

  //for -> while, eps 1e-12
  /*
    FCN=169.765 FROM MIGRAD    STATUS=CONVERGED     187 CALLS         188 TOTAL
                     EDM=8.97263e-10    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   1.0 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           1.99844e+05   4.51649e+02  -1.70227e-01  -1.19996e-07
   2  p1           1.00592e+02   3.60565e-02   2.20949e-05  -7.88267e-04
   3  p2           7.05192e+00   2.60869e-02  -2.64829e-06   6.48267e-06
   4  p3           9.89452e+00   4.04601e-02  -9.17655e-05  -7.11267e-04

   */

  //eps 1e-9
  /*
 FCN=169.765 FROM MIGRAD    STATUS=CONVERGED     187 CALLS         188 TOTAL
                     EDM=8.93042e-10    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   1.0 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  p0           1.99844e+05   4.51644e+02  -1.71038e-01  -1.19762e-07
   2  p1           1.00592e+02   3.60566e-02   2.20875e-05  -7.85793e-04
   3  p2           7.05192e+00   2.60868e-02  -2.70653e-06   7.86980e-06
   4  p3           9.89452e+00   4.04602e-02  -9.18923e-05  -7.07177e-04
   */
  const Double_t amp = TMath::Abs(rawamp);
  const Double_t mean = TMath::Abs(rawmean);
  const Double_t sigma = TMath::Abs(rawsigma);
  const Double_t lambda = TMath::Abs(rawlambda);

  const Double_t x0 = xx-mean;
  Double_t sum=0;
  Double_t sW = 0;

  //exp(-x/lambda)
  const Double_t step = 0.1;
  Double_t weight = 1;
  Double_t nlam = 0;
  const Double_t eps = 1e-9;
  while(weight>eps){
    weight = TMath::Exp(-nlam);
    const Double_t delta = TMath::Gaus( x0 - nlam*lambda, 0, sigma, kTRUE)*weight;
    sum += delta;
    sW += weight;
    nlam+=step;
  }
 
  if(sW>1e-12){
    return amp*sum/sW;
  }
  else{
    return 0;
  }
}

Double_t conv2(const Double_t xx, const Double_t amp, const Double_t mean, const Double_t sigma, const Double_t lambda)
{
  //lambda=0, exp(-x/lambda)/lambda -> delta function
  if(lambda<EPSILON){
    return amp*TMath::Gaus( xx, mean, sigma, kTRUE);
  }

  const Double_t kappa = 1/lambda;
  const Double_t eps = 1e-12;
  const Double_t nstep = 300;
  Double_t tau = 0;
  Double_t sum=0;
  for(Int_t ii=0; ii<nstep; ii++){
    sum+= TMath::Gaus(xx-tau, mean, sigma,kTRUE);
    
    //const Double_t oldx = tau;
    tau = -TMath::Log(TMath::Exp(-kappa*tau)-(1-eps/kappa)/nstep)/kappa;
    //printf("%d %e %e -- %e -- %e\n", ii,oldx, tau, exp(-kappa*tau)*kappa, exp(-kappa*oldx)*kappa-exp(-kappa*tau)*kappa);
  }
  
  return amp*sum/nstep;
}

/*
TF1 *f1=new TF1("f1","conv(x,0.1401767,49.556376,3.3698336, 52.608731)", 0,200)
f1->Draw()
 */

/*
int main(int argc, char * argv[])
{
  TList * ll = new TList;

  TH1D * hgau = new TH1D("hgau","",300,0,300); ll->Add(hgau);
  TH1D * hexp = new TH1D("hexp","",200,0,20); ll->Add(hexp);
  TH1D * htot = new TH1D("htot","",300,0,300); ll->Add(htot);

  const Int_t nsample = 2000000;

  const Double_t mean = 100;
  const Double_t sigma = 7;
  const Double_t lambda = 10.;

  TRandom3 ran(1);
  for(Int_t ii=0; ii<nsample; ii++){
    const Double_t xgau = ran.Gaus(mean, sigma);

    const Double_t GG = ran.Rndm();
    const Double_t xexp = -TMath::Log(1-GG)*lambda;

    hgau->Fill(xgau);
    hexp->Fill(xexp);
    htot->Fill(xgau+xexp);
  }

new TCanvas;
hh=htot
  hh->Draw();
gPad->SetLogy();
 
 TF1 *fg=new TF1("fg","gausn",0,300);
  TF1 *f1=new TF1("f1","conv(x,[0],[1],[2],[3])",0,300);
  hh->Fit(fg, "N");

  f1->SetParameters(fg->GetParameters());
  f1->SetParameter(3,fg->GetParameter(2)/10.);
  f1->SetLineColor(kBlue);  //f1->Draw("same");
  hh->Fit(f1, "N");

   fg->Draw("same");
  f1->Draw("same");

//===============
  new TCanvas;
  hgau->Draw();
  XGLUtils::SetFitPrintLevel(1);  
  XGLUtils::FFunc func = XGLUtils::TailedGaus; //XGLUtils::RawGaus; 
  Int_t npar = 4;
  Double_t pars[]={hgau->Integral(0,hgau->GetNbinsX())*hgau->GetBinWidth(1), hgau->GetMean(), hgau->GetRMS(),0};
  XGLUtils::ChisquareFit(hgau, func, npar, pars);  

  TH1D *hfitgau = XGLUtils::GetHfit(Form("%sfit", hgau->GetName()), func, pars, hgau->GetXaxis()->GetXmin(), hgau->GetXaxis()->GetXmax());
  hfitgau->SetLineColor(kRed);

  ll->Add(hfitgau);

  hfitgau->Draw("same");
//===



  XGLUtils::SetFitPrintLevel(1);  
  XGLUtils::FFunc func = XGLUtils::RawGaus; 
  Int_t npar = 3;
  Double_t pars[]={htot->Integral(0,htot->GetNbinsX())*htot->GetBinWidth(1), htot->GetMean(), htot->GetRMS(),0};
  XGLUtils::ChisquareFit(htot, func, npar, pars);  

  npar=4;
  func = XGLUtils::TailedGaus;
  pars[3]=1/(pars[2]/10);
  XGLUtils::ChisquareFit(htot, func, npar, pars);                   

  TH1D *hfit = XGLUtils::GetHfit(Form("%sfit", htot->GetName()), func, pars, htot->GetXaxis()->GetXmin(), htot->GetXaxis()->GetXmax());
  hfit->SetLineColor(kRed);

  ll->Add(hfit);

  hfit->Draw("same");

  TFile *fout=new TFile("outNonGaussian.root","recreate");
  ll->Write();
  fout->Save();
  fout->Close();

  return 0;
}
*/
